VIBRANT Trial - Primary Analyses

Author

Laura Symul, Laura Vermeren, Susan Holmes

Code

library(tidyverse)
library(magrittr)
library(gt)
library(patchwork)
library(SummarizedExperiment)
library(tidySummarizedExperiment)
library(brms)
library(gtsummary)
# library(mia) # BiocManager::install("mia")

theme_set(theme_light())
tmp <- fs::dir_map("R", source)

Code

data_source <- "simulated"

Data preparation

Loading the data and creating `SummarizedExperiment` objects for each assays

In this section, we load the raw data from their original files (i.e., the data/files as generated by the corresponding labs or teams) and transform each assay into a SummarizedExperiment object. This will enable to then collect them into a single MultiAssayExperiment object and perform the integrative analyses.

To link the assays together,we define a common unique identifier for each participant x visit. This unique identifier (.sample) is the concatenation of the pid (participant ID) and the visit_code.

Clinical and survey (CRF) data

Information about the data

The data acquisition process is described XXXHEREXXXX.

The data encoding in a xxx database and quality control processes are described XXXHEREXXX.

The data was exported from the database by Lara Lewis (CAPRISA, South Africa) as a set of xxx files.

Lara Wautier (UCLouvain, Belgium) then imported these files in R and assembled the data into two tables: a participants and a visits table. The code is available XXXHEREXXX.

The participants table contains participant-invariant information such as the participant ID, the treatment arm, site, and study population (ITT, mITT, PP). The visits table contains participant-varying information collected at each visit.

The columns necessary for the primary trial analyses were then selected from these two and exported into two files: participants.RDS and visits.RDS, which we load below.

Loading the data

We load the participant and visits data, and create

one table which will become the MAE’s colData
one SE object which provides attendance information at each visit (this also ensures that the MAE object has data for all visits, which is useful for downstream analyses)

Code

participants <- readRDS(file = str_c(data_dir(), "00 clinical data/participants.RDS"))
visits <- readRDS(file = str_c(data_dir(), "00 clinical data/visits.RDS"))

clin <- 
  participants |> left_join(visits, by = join_by(pid)) |> 
  mutate(sample_id = str_c(pid, "_", visit_code), rownames = sample_id) |> 
  dplyr::select(sample_id, everything()) 

mae_coldata <- 
  clin |> 
  as.data.frame() |> 
  column_to_rownames(var = "rownames")

Code

make_se_clin <- function(clin, cols, assayname){
  assay <- 
    clin |> 
    dplyr::select(sample_id, all_of(cols)) |> 
    as.data.frame() |> 
    column_to_rownames("sample_id") |> 
    t()
  
  assay_list <- list()
  assay_list[[assayname]] <- assay
  
  SummarizedExperiment::SummarizedExperiment(
    assays = assay_list
  )
  
}

Code

# Directly included in data simulation 

# clin <- clin |> mutate(planned = TRUE, attended = ifelse(str_detect(visit_name, "Daily"), NA, attended)) #TRUE

Code

se_va <- make_se_clin(clin, cols = c("planned", "attended"), assayname = "visit_attendance")

Weekly metagenomics data with LBP strains relative abundances

TODO: add information about how the data was generated (summary of the methods)

We load the weekly metagenomics data, the manifest for the metagenomic barcodes, and the strain information.

Code

mg <- read_csv(str_c(data_dir(), "01 metagenomics kSanity/mg_kSanity.csv"))
mg_combined <- read_csv(str_c(data_dir(), "03 metagenomics combined/mg_combined.csv"))
mg_manifest <- read_csv(str_c(data_dir(), "01 metagenomics kSanity/mg_kSanity_manifest.csv"))
LBP_strain_info <- readxl::read_xlsx(str_c(dropbox_dir(), "9_VIBRANT_Study Product/IsolateNumbers.xlsx"))

LBP_strain_info <- 
 LBP_strain_info |> 
  mutate(strain_id = `VMRC ID`, 
         LBP = ifelse(is.na(LC106), "LC-115", "LC-106 & LC-115") |> factor(), 
         strain_origin = `Geographic source` |> factor()
  ) |> 
  dplyr::rename(Biose_ID = `Biose ID`, VMRC_ID = `VMRC ID`) |> 
  arrange(strain_origin, LBP) |> 
  mutate(strain_id = strain_id |> fct_inorder()) |> 
  dplyr::select(strain_id, LBP, strain_origin, contains("ID"))

LBP_strain_info |> gt(caption = "LBP strain information")

LBP strain information
strain_id	LBP	strain_origin	Biose_ID	VMRC_ID
FF00018	LC-106 & LC-115	SA	GA08	FF00018
FF00051	LC-106 & LC-115	SA	GA09	FF00051
UC101_2_33	LC-106 & LC-115	SA	GA12	UC101_2_33
FF00004	LC-115	SA	GA07	FF00004
FF00064	LC-115	SA	GA10	FF00064
FF00072	LC-115	SA	GA11	FF00072
UC119_2_PB0238	LC-115	SA	GA13	UC119_2_PB0238
122010_1999_16	LC-115	SA	GA14	122010_1999_16
185329_1999_17	LC-115	SA	GA15	185329_1999_17
C0059E1	LC-106 & LC-115	US	GA03	C0059E1
C0022A1	LC-106 & LC-115	US	GA04	C0022A1
C0175A1	LC-106 & LC-115	US	GA06	C0175A1
C0006A1	LC-115	US	GA01	C0006A1
CC0028A1	LC-115	US	GA02	CC0028A1
C0112A1	LC-115	US	GA05	C0112A1

Code

mg_to_SE <- function(mg,mg_combined, mg_manifest, LBP_strain_info){
  mg_manifest <- 
    mg_manifest |> 
    mutate(sample_id = str_c(pid, "_", visit_code)) # nécéssaire ? pourquoi pas "barcode" ccréer dans sim ? 
  
  if (any(duplicated(mg_manifest$sample_id))) stop("Duplicated `uid` in `mg_manifest`")
  
  # mg <- mg |> filter(!is.na(`Lactobacillus crispatus`))
  mg <- mg |> left_join(mg_manifest, by = join_by(barcode))
  # mg_combined <- mg_combined |> filter(!is.na(`Lactobacillus crispatus`))
  mg_combined <- mg_combined |> left_join(mg_manifest, by = join_by(barcode))
  tmp <- mg |> dplyr::select(barcode, sample_id) |> distinct()
  if (any(duplicated(tmp$sample_id))) stop("Duplicated `uid` in `mg`")
  
  assay_prop_of_Lc <- 
    mg |> 
    # dplyr::select(sample_id, strain, prop_of_Lc) |> 
    # pivot_wider(names_from = sample_id, values_from = prop_of_Lc, values_fill = 0) |>  # the new database is already in wide format 
    dplyr::select(sample_id, everything()) |> 
    dplyr::select(-c(pid, visit_code, barcode)) |> 
    as.data.frame() |> 
    column_to_rownames("sample_id") |> 
    t() 
  # colnames(assay_prop_of_Lc) <- assay_prop_of_Lc[1, ]  
  # assay_prop_of_Lc <- assay_prop_of_Lc[-1, ]  
    
  assay_unique_kmers <- 
    mg_combined |> 
    # dplyr::select(sample_id, strain, unique_kmers) |> 
    # pivot_wider(names_from = sample_id, values_from = unique_kmers, values_fill = 0) |> 
    dplyr::select(sample_id, all_of(LBP_strain_info$strain_id), `Lactobacillus crispatus`) |> 
    as.data.frame() |>
    column_to_rownames("sample_id") |>
    t()

  se_coldata <-
    mg |>
    dplyr::select(sample_id, barcode) |> #shared_kmers
    distinct() |>
    as.data.frame()

  se_rowdata <- 
    mg |> 
    dplyr::select(-c(barcode, pid, visit_code, sample_id)) |> 
    colnames() |> 
    as.data.frame() |> 
    setNames("strain") |> 
    distinct() |> 
    left_join(
      LBP_strain_info,
      by = join_by(strain == strain_id)
    ) 
  
  
  # Harmonization of the order of samples and feature
  
  sorted_sample_ids <- sort(as.character(se_coldata$sample_id))
  assay_prop_of_Lc <- assay_prop_of_Lc[, sorted_sample_ids]
  assay_unique_kmers <- assay_unique_kmers[, sorted_sample_ids]
  se_coldata <- se_coldata[match(sorted_sample_ids, se_coldata$sample_id), ]

  sorted_strains <- assay_prop_of_Lc |> rownames() # sort(as.character(se_rowdata$strain))
  assay_prop_of_Lc <- assay_prop_of_Lc[sorted_strains, ]
  assay_unique_kmers <- assay_unique_kmers[sorted_strains, ]
  se_rowdata <- se_rowdata[match(sorted_strains, se_rowdata$strain), ]

  se_rowdata <- as(se_rowdata, "DataFrame")
  se_coldata <- as(se_coldata, "DataFrame")

  SummarizedExperiment::SummarizedExperiment(
    assays = list(prop_of_Lc = assay_prop_of_Lc, unique_kmers = assay_unique_kmers),
    rowData = se_rowdata,
    colData = se_coldata
  )
}

Code

(se_mg <- mg_to_SE(mg,mg_combined,  mg_manifest, LBP_strain_info))

# A SummarizedExperiment-tibble abstraction: 14,832 × 11
# Features=16 | Samples=927 | Assays=prop_of_Lc, unique_kmers
   .feature       .sample prop_of_Lc unique_kmers sample_id barcode strain LBP  
   <chr>          <chr>        <dbl>        <dbl> <chr>     <chr>   <chr>  <fct>
 1 FF00018        p1_1000          0            0 p1_1000   mg_p1_… FF000… LC-1…
 2 FF00051        p1_1000          0            0 p1_1000   mg_p1_… FF000… LC-1…
 3 UC101_2_33     p1_1000          0            0 p1_1000   mg_p1_… UC101… LC-1…
 4 C0059E1        p1_1000          0            0 p1_1000   mg_p1_… C0059… LC-1…
 5 C0022A1        p1_1000          0            0 p1_1000   mg_p1_… C0022… LC-1…
 6 C0175A1        p1_1000          0            0 p1_1000   mg_p1_… C0175… LC-1…
 7 FF00004        p1_1000          0            0 p1_1000   mg_p1_… FF000… LC-1…
 8 FF00064        p1_1000          0            0 p1_1000   mg_p1_… FF000… LC-1…
 9 FF00072        p1_1000          0            0 p1_1000   mg_p1_… FF000… LC-1…
10 UC119_2_PB0238 p1_1000          0            0 p1_1000   mg_p1_… UC119… LC-1…
# ℹ 310 more rows
# ℹ 3 more variables: strain_origin <fct>, Biose_ID <chr>, VMRC_ID <chr>

Checks & Quality Controls

Total unique kmers per sample

Code

se_mg |> 
  as_tibble() |> 
  mutate(unique_kmers = as.numeric(unique_kmers)) |>
  group_by(.sample) |> 
  summarize(total_unique_kmers = sum(unique_kmers)) |> 
  ggplot(aes(x = total_unique_kmers)) +
  geom_histogram() +
  scale_x_log10()

Data overview

Code

se_mg |> 
  as_tibble() |> 
  mutate(.feature = .feature |> fct_inorder(),
         prop_of_Lc = as.numeric(prop_of_Lc)) |> 
  ggplot(aes(x = .feature, y = .sample, fill = prop_of_Lc)) +
  geom_tile() +
  scale_fill_continuous(low = "white", high = "steelblue2") +
  scale_x_discrete("Strains") +
  scale_y_discrete("Samples",breaks = NULL) +
  facet_grid(. ~ LBP + strain_origin, scales = "free_x", space = "free") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Weekly 16S rRNA data

TODO: add a summary of the protocols and methods used to generate the 16S data

We load the 16S rRNA data, which is stored as a phyloseq object.

Code

(w16S <- readRDS(str_c(data_dir(), "03 weekly 16S/w16S.RDS")))

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 20 taxa and 927 samples ]
sample_data() Sample Data:       [ 927 samples by 2 sample variables ]
tax_table()   Taxonomy Table:    [ 20 taxa by 4 taxonomic ranks ]

We transform the phyloseq object into a SummarizedExperiment object.

Code

phyloseq_to_SE <- function(physeq){
  
  sample_info <- 
    physeq@sam_data |>
    as.data.frame() |> 
    set_class("data.frame") |> 
    rownames_to_column("barcode") |> 
    as_tibble() |> 
    mutate(sample_id = str_c(pid, "_", visit_code))

  counts <- physeq@otu_table |> as.data.frame() |> set_class("data.frame")
  if (!phyloseq::taxa_are_rows(physeq)) counts <- counts |> t()
  counts <- counts |> set_colnames(sample_info$sample_id)
  rel <- t(t(counts)/colSums(counts)) |> as.data.frame()
  
  SummarizedExperiment::SummarizedExperiment(
    assays = list(ASV_counts = counts, ASV_prop = rel),
    rowData = physeq@tax_table |> as.data.frame() |> set_class("data.frame"),
    colData = sample_info |> as.data.frame() |> column_to_rownames("sample_id")
  )
}

Code

(se_16S <- phyloseq_to_SE(w16S))

# A SummarizedExperiment-tibble abstraction: 18,540 × 11
# Features=20 | Samples=927 | Assays=ASV_counts, ASV_prop
   .feature   .sample ASV_counts ASV_prop barcode pid   visit_code asv   Species
   <chr>      <chr>        <dbl>    <dbl> <chr>   <fct>      <dbl> <chr> <chr>  
 1 Lactobaci… p1_1000          0 0        16S_p1… p1          1000 "Lac… iners  
 2 Lactobaci… p1_1000         37 0.0300   16S_p1… p1          1000 "Lac… crispa…
 3 Lactobaci… p1_1000          0 0        16S_p1… p1          1000 "Lac… gasseri
 4 Lactobaci… p1_1000          0 0        16S_p1… p1          1000 "Lac… jensen…
 5 Gardnerel… p1_1000         18 0.0146   16S_p1… p1          1000 "Gar… leopol…
 6 Gardnerel… p1_1000          0 0        16S_p1… p1          1000 "Gar… swidsi…
 7 Gardnerel… p1_1000        260 0.211    16S_p1… p1          1000 "Gar… vagina…
 8 Bifidobac… p1_1000         77 0.0623   16S_p1… p1          1000 "Bif… breve  
 9 Bifidobac… p1_1000        112 0.0907   16S_p1… p1          1000 "Bif… longum 
10 Megasphae… p1_1000          1 0.000810 16S_p1… p1          1000 "Meg… lornae 
# ℹ 390 more rows
# ℹ 2 more variables: is_lacto <chr>, color <chr>

Checks & Quality Controls

16S Total reads per sample

Code

se_16S |> 
  as_tibble() |> 
  group_by(.sample) |> summarize(total_counts = sum(ASV_counts)) |> 
  ggplot(aes(x = total_counts)) +
  geom_histogram() +
  scale_x_log10()

Data overview

Code

se_16S |> 
  as_tibble() |> 
  group_by(.feature) |> mutate(total_prop = sum(ASV_prop)) |> ungroup() |> 
  arrange(-total_prop) |> mutate(.feature = .feature |> fct_inorder()) |> 
  ggplot(aes(x = .feature, y = .sample, fill = asinh(ASV_counts))) +
  geom_tile() +
  scale_fill_continuous(low = "white", high = "steelblue2") +
  scale_x_discrete("ASVs", breaks = NULL) +
  scale_y_discrete("Samples",breaks = NULL)

Integrating all assays into a `MultiAssayExperiment` object

In this section, we assemble all SE objects created in the previous section into a single MultiAssayExperiment object.

`colData`

Code

mae_coldata |> as_tibble()

# A tibble: 4,738 × 20
   sample_id site  arm   pid     age race  food_insecurity contraception      
   <chr>     <fct> <fct> <fct> <int> <chr> <chr>           <chr>              
 1 p1_0      CAP   Pl    p1       25 White No              Non-hormonal method
 2 p1_1000   CAP   Pl    p1       25 White No              Non-hormonal method
 3 p1_1001   CAP   Pl    p1       25 White No              Non-hormonal method
 4 p1_1002   CAP   Pl    p1       25 White No              Non-hormonal method
 5 p1_1003   CAP   Pl    p1       25 White No              Non-hormonal method
 6 p1_1004   CAP   Pl    p1       25 White No              Non-hormonal method
 7 p1_1005   CAP   Pl    p1       25 White No              Non-hormonal method
 8 p1_1006   CAP   Pl    p1       25 White No              Non-hormonal method
 9 p1_1007   CAP   Pl    p1       25 White No              Non-hormonal method
10 p1_1100   CAP   Pl    p1       25 White No              Non-hormonal method
# ℹ 4,728 more rows
# ℹ 12 more variables: gender_sexual_partners <chr>,
#   nb_partner_pastmonth <dbl>, nb_partner_life <dbl>, dropout <lgl>,
#   dropout_time <int>, visit_code <dbl>, visit_abbr <fct>, visit_name <fct>,
#   visit_number <dbl>, day <dbl>, attended <lgl>, planned <lgl>

Assays

Visit attendance

Code

se_va

# A SummarizedExperiment-tibble abstraction: 9,476 × 3
# Features=2 | Samples=4738 | Assays=visit_attendance
   .feature .sample visit_attendance
   <chr>    <chr>   <lgl>           
 1 planned  p1_0    TRUE            
 2 attended p1_0    FALSE           
 3 planned  p1_1000 TRUE            
 4 attended p1_1000 TRUE            
 5 planned  p1_1001 TRUE            
 6 attended p1_1001 NA              
 7 planned  p1_1002 TRUE            
 8 attended p1_1002 NA              
 9 planned  p1_1003 TRUE            
10 attended p1_1003 NA              
# ℹ 30 more rows

Weekly metagenomics data with LBP strains relative abundances

Code

se_mg

# A SummarizedExperiment-tibble abstraction: 14,832 × 11
# Features=16 | Samples=927 | Assays=prop_of_Lc, unique_kmers
   .feature       .sample prop_of_Lc unique_kmers sample_id barcode strain LBP  
   <chr>          <chr>        <dbl>        <dbl> <chr>     <chr>   <chr>  <fct>
 1 FF00018        p1_1000          0            0 p1_1000   mg_p1_… FF000… LC-1…
 2 FF00051        p1_1000          0            0 p1_1000   mg_p1_… FF000… LC-1…
 3 UC101_2_33     p1_1000          0            0 p1_1000   mg_p1_… UC101… LC-1…
 4 C0059E1        p1_1000          0            0 p1_1000   mg_p1_… C0059… LC-1…
 5 C0022A1        p1_1000          0            0 p1_1000   mg_p1_… C0022… LC-1…
 6 C0175A1        p1_1000          0            0 p1_1000   mg_p1_… C0175… LC-1…
 7 FF00004        p1_1000          0            0 p1_1000   mg_p1_… FF000… LC-1…
 8 FF00064        p1_1000          0            0 p1_1000   mg_p1_… FF000… LC-1…
 9 FF00072        p1_1000          0            0 p1_1000   mg_p1_… FF000… LC-1…
10 UC119_2_PB0238 p1_1000          0            0 p1_1000   mg_p1_… UC119… LC-1…
# ℹ 310 more rows
# ℹ 3 more variables: strain_origin <fct>, Biose_ID <chr>, VMRC_ID <chr>

Weekly 16S rRNA data

Code

se_16S

# A SummarizedExperiment-tibble abstraction: 18,540 × 11
# Features=20 | Samples=927 | Assays=ASV_counts, ASV_prop
   .feature   .sample ASV_counts ASV_prop barcode pid   visit_code asv   Species
   <chr>      <chr>        <dbl>    <dbl> <chr>   <fct>      <dbl> <chr> <chr>  
 1 Lactobaci… p1_1000          0 0        16S_p1… p1          1000 "Lac… iners  
 2 Lactobaci… p1_1000         37 0.0300   16S_p1… p1          1000 "Lac… crispa…
 3 Lactobaci… p1_1000          0 0        16S_p1… p1          1000 "Lac… gasseri
 4 Lactobaci… p1_1000          0 0        16S_p1… p1          1000 "Lac… jensen…
 5 Gardnerel… p1_1000         18 0.0146   16S_p1… p1          1000 "Gar… leopol…
 6 Gardnerel… p1_1000          0 0        16S_p1… p1          1000 "Gar… swidsi…
 7 Gardnerel… p1_1000        260 0.211    16S_p1… p1          1000 "Gar… vagina…
 8 Bifidobac… p1_1000         77 0.0623   16S_p1… p1          1000 "Bif… breve  
 9 Bifidobac… p1_1000        112 0.0907   16S_p1… p1          1000 "Bif… longum 
10 Megasphae… p1_1000          1 0.000810 16S_p1… p1          1000 "Meg… lornae 
# ℹ 390 more rows
# ℹ 2 more variables: is_lacto <chr>, color <chr>

Assembling the `MultiAssayExperiment` object

Code

assay_list <- 
  list(
    visit_attendance = se_va,
    mg_ksanity = se_mg,
    # mg_virgo2 = se_mg_V,
    amplicon_ASV_all = se_16S
    )

mae <- 
  MultiAssayExperiment::MultiAssayExperiment(
    experiments = MultiAssayExperiment::ExperimentList(assay_list),
    colData = mae_coldata |> dplyr::rename(.sample = sample_id),
    metadata = list(data_source = data_source, date_created = today())
  )

Code

mae

A MultiAssayExperiment object of 3 listed
 experiments with user-defined names and respective classes.
 Containing an ExperimentList class object of length 3:
 [1] visit_attendance: SummarizedExperiment with 2 rows and 4738 columns
 [2] mg_ksanity: SummarizedExperiment with 16 rows and 927 columns
 [3] amplicon_ASV_all: SummarizedExperiment with 20 rows and 927 columns
Functionality:
 experiments() - obtain the ExperimentList instance
 colData() - the primary/phenotype DataFrame
 sampleMap() - the sample coordination DataFrame
 `$`, `[`, `[[` - extract colData columns, subset, or experiment
 *Format() - convert into a long or wide DataFrame
 assays() - convert ExperimentList to a SimpleList of matrices
 exportClass() - save data to flat files

Checks, QC, and Outcomes Definitions

Checks and Quality Controls

Visualization of available data for each assay and each `site x arm x participant x visit`

Code

available_data <- 
  mae[["visit_attendance"]] |> 
  as_tibble() |> 
  filter(.feature == "attended") |> 
  dplyr::rename(attended = visit_attendance) |> 
  left_join(
    mae@sampleMap |> 
      as_tibble() |> 
      filter(assay != "visit_attendance") |> 
      mutate(.sample = primary) |> 
      dplyr::select(.sample, assay),
    by = join_by(.sample)
  ) |> 
  left_join(
    mae@colData |> as_tibble() |> dplyr::select(.sample, pid, site, arm, starts_with("visit_"))
  ) |> 
  mutate(available_data = ifelse(attended == FALSE, FALSE, TRUE))

assay_collection_schedule <- 
  available_data |> 
  filter(!is.na(assay)) |> 
  group_by(assay, visit_code) |> 
  summarize(assay_scheduled_at_visit = !all(is.na(available_data)), .groups = "drop") |> 
  filter(assay_scheduled_at_visit)

available_data <- 
  available_data |> 
  left_join(
    assay_collection_schedule,
    by = join_by(assay, visit_code)
  )

available_data <- 
  available_data |> 
  filter(assay_scheduled_at_visit) 

available_data

# A tibble: 1,854 × 13
   .feature .sample attended assay       pid   site  arm   visit_code visit_abbr
   <chr>    <chr>   <lgl>    <fct>       <fct> <fct> <fct>      <dbl> <fct>     
 1 attended p1_1000 TRUE     mg_ksanity  p1    CAP   Pl          1000 ENR       
 2 attended p1_1000 TRUE     amplicon_A… p1    CAP   Pl          1000 ENR       
 3 attended p1_1100 FALSE    mg_ksanity  p1    CAP   Pl          1100 Post-ABX  
 4 attended p1_1100 FALSE    amplicon_A… p1    CAP   Pl          1100 Post-ABX  
 5 attended p1_1200 TRUE     mg_ksanity  p1    CAP   Pl          1200 WFU       
 6 attended p1_1200 TRUE     amplicon_A… p1    CAP   Pl          1200 WFU       
 7 attended p1_1300 TRUE     mg_ksanity  p1    CAP   Pl          1300 WFU       
 8 attended p1_1300 TRUE     amplicon_A… p1    CAP   Pl          1300 WFU       
 9 attended p1_1400 TRUE     mg_ksanity  p1    CAP   Pl          1400 WFU       
10 attended p1_1400 TRUE     amplicon_A… p1    CAP   Pl          1400 WFU       
# ℹ 1,844 more rows
# ℹ 4 more variables: visit_name <fct>, visit_number <dbl>,
#   available_data <lgl>, assay_scheduled_at_visit <lgl>

Code

available_data |> 
  filter(attended == TRUE) |>  # ou available_data == TRUE
  ggplot(aes(x = assay, y = pid, col = assay)) +
  geom_point() +
  facet_grid(arm + site ~ visit_code, scales = "free", space = "free") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), legend.position  = "none")

Checking for contamination

While it is possible for a few Placebo participants to natively have strains that are extremely similar to the LBP strains, it is highly unlikely for many of the LBP strains to be present simultaneously in a single sample for a placebo participant. Such event would much more likely be the result of contamination OR a mislabeling of the sample OR a randomization issue.

We thus check how many LBP strains are simultaneously present in the Placebo participants samples

TODO

Code

mae[["mg_ksanity"]] |> 
  left_join(
    colData(mae) |> 
      as.tibble() |> 
      dplyr::select(.sample, site, arm, visit_code, pid), 
    by = ".sample"
  ) |> 
  filter(arm == "Pl" & !is.na(LBP)) |> 
  ggplot(aes(x = pid, y = unique_kmers, color = strain)) + 
  geom_point() + 
  facet_grid(visit_code ~ strain, scales = "free", space = "free") + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))

Other checks ?

Outcomes definitions

Primary outcome definition

The primary outcome is XXX INSERT DEFINITION XXX

Consequently, we need to compute the relative abundance of each LBP strains from their proportions out of total L. crispatus proportions and the total L. crispatus proportions as estimated by 16S sequencing.

LBP strains relative abundances in each samples

This is stored in a new assay: LBP_rel_ab

Code

se_LBP_rel_abundances <- 
  list(
    LBP_rel_ab = 
      mae[["mg_ksanity"]] |> 
      filter(!is.na(LBP)) |> 
      left_join(
        mae[["amplicon_ASV_all"]] |> 
          filter(str_detect(.feature, "crispatus")) |> 
          as_tibble() |> 
          dplyr::select(.sample, ASV_prop) |> 
          group_by(.sample) |> summarize(Lc_prop = sum(ASV_prop))
      ) |> 
      mutate(
        LBP_rel_ab = prop_of_Lc * Lc_prop
      ) |> 
      dplyr::select(-prop_of_Lc, -unique_kmers) 
  )


mae <- c(mae, se_LBP_rel_abundances)

At each visit

From this, we can compute the total relative abundance of all LBP strains or the maximum relative abundance of any LBP strain in each sample.

If the total relative abundance of all LBP strains is larger than 0.1 or the maximum relative abundance of any LBP strain is larger than 0.05, we consider that the LBP achieved colonization at that visit.

Code

colonization <- 
  mae[["LBP_rel_ab"]] |> 
  as_tibble() |>
  left_join(colData(mae) |> as_tibble(),by = join_by(.sample)) |> 
  group_by(.sample, pid, visit_code) |> 
  summarize(
    total_prop_at = sum(LBP_rel_ab), # , na.rm = TRUE
    max_prop_at = max(LBP_rel_ab), # , na.rm = TRUE
    .groups = "drop"
  ) |> 
  mutate(
    colonized_at = (total_prop_at > 0.1) | max_prop_at > 0.05
  )

By each visit

From the colonization status at each visit, we can compute the colonization status by each visit. A participant is considered to have colonized by a visit if they have been colonized at that visit or any previous visit.

Code

colonization <- 
  colonization |> 
  group_by(pid) |> 
  arrange(pid, visit_code) |>
  mutate(
    colonized_by = cummax(ifelse(is.na(colonized_at), FALSE, colonized_at)) |> as.logical()
  )  |> 
  
  # mutate(
  #   colonized_by = cumsum(colonized_at) > 0
  # ) |> 
  
  ungroup()

Adding `primary` assay to `MAE`

Code

se_primary <-
  SummarizedExperiment(
    assays = list(
      primary = 
        colonization |> 
        dplyr::select(-pid, -visit_code) |> 
        as.data.frame() |> 
        column_to_rownames(".sample") |> t()
    )
  )

mae <- c(mae, list(LBP_colonization = se_primary))

Secondary outcomes definitions

16S rRNA microbiota composition, aggregated at the species level

TODO

16S rRNA microbiota composition, augmented with LBP strains relative abundances

We combine the 16S rRNA microbiota composition with the LBP strains relative abundances to create a new assay mb_composition.

TODO

We now have an “augmented” MAE object with the following assays:

Code

mae

A MultiAssayExperiment object of 5 listed
 experiments with user-defined names and respective classes.
 Containing an ExperimentList class object of length 5:
 [1] visit_attendance: SummarizedExperiment with 2 rows and 4738 columns
 [2] mg_ksanity: SummarizedExperiment with 16 rows and 927 columns
 [3] amplicon_ASV_all: SummarizedExperiment with 20 rows and 927 columns
 [4] LBP_rel_ab: SummarizedExperiment with 15 rows and 927 columns
 [5] LBP_colonization: SummarizedExperiment with 4 rows and 927 columns
Functionality:
 experiments() - obtain the ExperimentList instance
 colData() - the primary/phenotype DataFrame
 sampleMap() - the sample coordination DataFrame
 `$`, `[`, `[[` - extract colData columns, subset, or experiment
 *Format() - convert into a long or wide DataFrame
 assays() - convert ExperimentList to a SimpleList of matrices
 exportClass() - save data to flat files

Analyses

Demographics (Table 1)

In this section we will create a table with the demographics of the participants in the study.

Code

participants |> 
  select(-c("pid", "dropout_time"))  |> 
  tbl_summary(
    by = arm, 
    type = 
      list(
        nb_partner_pastmonth ~ "continuous", 
        nb_partner_life ~ "continuous"
      )
  )  |>
  add_p(include = -c("site"))

Characteristic	Pl N = 23¹	LC-106-7 N = 27¹	LC-106-3 N = 15¹	LC-106-o N = 16¹	LC-115 N = 22¹	p-value²
site
CAP	11 (48%)	15 (56%)	11 (73%)	12 (75%)	11 (50%)
MGH	12 (52%)	12 (44%)	4 (27%)	4 (25%)	11 (50%)
age	26 (22, 33)	27 (24, 31)	26 (19, 34)	26 (20, 34)	29 (21, 36)	>0.9
race
Asian	3 (13%)	2 (7.4%)	5 (33%)	2 (13%)	5 (23%)
Black	9 (39%)	6 (22%)	0 (0%)	2 (13%)	7 (32%)
Multiple	8 (35%)	9 (33%)	6 (40%)	9 (56%)	6 (27%)
White	3 (13%)	10 (37%)	4 (27%)	3 (19%)	4 (18%)
food_insecurity	5 (22%)	10 (37%)	4 (27%)	6 (38%)	8 (36%)	0.7
contraception
Combined oral contraceptive pills	4 (17%)	4 (15%)	1 (6.7%)	0 (0%)	2 (9.1%)
Copper IUD	2 (8.7%)	9 (33%)	2 (13%)	3 (19%)	4 (18%)
Depo Provera	3 (13%)	1 (3.7%)	2 (13%)	2 (13%)	3 (14%)
Hormone containing IUD	5 (22%)	5 (19%)	6 (40%)	3 (19%)	4 (18%)
Net EN	2 (8.7%)	4 (15%)	2 (13%)	7 (44%)	7 (32%)
Non-hormonal method	7 (30%)	4 (15%)	2 (13%)	1 (6.3%)	2 (9.1%)
gender_sexual_partners						0.4
Both	7 (30%)	7 (26%)	6 (40%)	10 (63%)	7 (32%)
Female	11 (48%)	11 (41%)	6 (40%)	3 (19%)	7 (32%)
Male	5 (22%)	9 (33%)	3 (20%)	3 (19%)	8 (36%)
nb_partner_pastmonth	1.00 (0.00, 1.00)	0.00 (0.00, 1.00)	0.00 (0.00, 1.00)	0.00 (0.00, 0.50)	0.50 (0.00, 2.00)	0.2
nb_partner_life	3.0 (2.0, 5.0)	3.0 (2.0, 5.0)	2.0 (2.0, 4.0)	2.5 (1.5, 4.0)	2.5 (0.0, 6.0)	0.6
dropout	3 (13%)	7 (26%)	1 (6.7%)	2 (13%)	2 (9.1%)	0.5
¹ n (%); Median (Q1, Q3)
² Kruskal-Wallis rank sum test; Fisher’s exact test

Primary outcomes (Table 2)

Colonization by week 5

Code

#simulation
col_by_week5 <- 
  expand_grid(
    site = c("MGH", "CAP") |> factor(),
    arm = c("Pl", "LC-106-7", "LC-106-3", "LC-106-o", "LC-115") |> fct_inorder(),
    pid = 1:10
  ) |> 
  slice_sample(prop = 0.9) |> 
  mutate(
    p = 0.0 + c(0, 0.5, 0.4, 0.7, 0.6)[arm |> as.integer()],
    colonized_by = (rbinom(n(), 1, p) == 1)
  )

col_by_week5 <- 
  col_by_week5 |> 
  mutate(out = (site == "MGH") & (arm %in% c("LC-106-3", "LC-106-o")) & (pid > 4)) |>
  filter(!out)

Code

col_by_week5 <-
  colonization |>
  left_join(colData(mae) |> as_tibble()) |>
  filter(visit_code == 1500)

Visualization of the primary outcome

Code

col_by_week5 |> 
  arrange(site, arm, colonized_by) |>
  group_by(site, arm) |> 
  mutate(participant_number = row_number() |> factor()) |> 
  ungroup() |> 
  ggplot(aes(x = participant_number, y = site, col = site, fill = site, shape = colonized_by)) +
  geom_point(size = 3) +
  facet_grid(. ~ arm, scales = "free_x") +
  guides(color = "none", fill = "none") +
  scale_shape_manual("Colonized by LBP by week 5", values = c(1, 16)) +
  xlab("Participant number") + ylab("")

Code

summary_table <- 
  col_by_week5 |> 
  group_by(site, arm) |> 
  summarize(
    n = n(), 
    n_success = sum(colonized_by),
    .groups = "drop"
  ) |> 
  mutate(
    p = n_success / n,
    CI = binom::binom.confint(n_success, n, method = "exact")
  )

Code

summary_table |> 
  mutate(
    N_per_site_and_arm = str_c("N = ", n),
    LBP_strain_detected = str_c(n_success, " (", round(p * 100)," %)"),
    CI_text = str_c(round(CI$lower * 100), "% - ", round(CI$upper * 100), "%")
  ) |> 
  dplyr::select(
    site, arm, N_per_site_and_arm, LBP_strain_detected, CI_text
  ) |> 
  pivot_longer(-c(site, arm)) |>
  mutate(
    name = 
      case_when(
        name == "N_per_site_and_arm" ~ "N participants",
        name == "LBP_strain_detected" ~ "n (%) of participants\nwith LBP strain detected",
        name == "CI_text" ~ "95% CI"
      )
  ) |> 
  pivot_wider(id_cols = c(site, name), names_from = arm, values_from = value) |> 
  group_by(site) |> 
  gt(caption = "Table 1: Colonization by week 5 by arm and site", row_group_as_column = TRUE) |> 
  cols_width(
    "name" ~ px(200),
    everything() ~ px(120)
  ) |> 
  cols_label(
    name = "",
    Pl = "Placebo",
    `LC-106-7` = "LC-106<br>7 days",
    `LC-106-3` = "LC-106<br>3 days",
    `LC-106-o` = "LC-106<br>overlap",
    `LC-115` = "LC-115<br>7 days",
    .fn = md
  )

Table 1: Colonization by week 5 by arm and site
		Placebo	LC-106 7 days	LC-106 3 days	LC-106 overlap	LC-115 7 days
CAP	N participants	N = 11	N = 15	N = 11	N = 12	N = 11
	n (%) of participants with LBP strain detected	0 (0 %)	5 (33 %)	8 (73 %)	11 (92 %)	7 (64 %)
	95% CI	0% - 28%	12% - 62%	39% - 94%	62% - 100%	31% - 89%
MGH	N participants	N = 12	N = 12	N = 4	N = 4	N = 11
	n (%) of participants with LBP strain detected	0 (0 %)	8 (67 %)	3 (75 %)	4 (100 %)	5 (45 %)
	95% CI	0% - 26%	35% - 90%	19% - 99%	40% - 100%	17% - 77%

TODO: include the p-values from model/test below

Tests / models

If we have 0 successes in none of the Placebo arm, we won’t be able to use a logistic regression model that accounts for site differences since we have a “perfect separation issue” and the odds are infinite for the Placebo arm:

Code

fit <- glm(colonized_by ~ arm + site + arm:site, data = col_by_week5, family = binomial)
fit |> summary()


Call:
glm(formula = colonized_by ~ arm + site + arm:site, family = binomial, 
    data = col_by_week5)

Coefficients:
                      Estimate Std. Error z value Pr(>|z|)
(Intercept)         -1.857e+01  1.967e+03  -0.009    0.992
armLC-106-7          1.787e+01  1.967e+03   0.009    0.993
armLC-106-3          1.955e+01  1.967e+03   0.010    0.992
armLC-106-o          2.096e+01  1.967e+03   0.011    0.991
armLC-115            1.913e+01  1.967e+03   0.010    0.992
siteMGH              7.888e-11  2.723e+03   0.000    1.000
armLC-106-7:siteMGH  1.386e+00  2.723e+03   0.001    1.000
armLC-106-3:siteMGH  1.178e-01  2.723e+03   0.000    1.000
armLC-106-o:siteMGH  1.617e+01  4.248e+03   0.004    0.997
armLC-115:siteMGH   -7.419e-01  2.723e+03   0.000    1.000

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 142.779  on 102  degrees of freedom
Residual deviance:  88.224  on  93  degrees of freedom
AIC: 108.22

Number of Fisher Scoring iterations: 17

There are a few ways to remedy that:

One way is to use a “pseudo-Bayesian” approach and to add 1 artificial success to all arms, but this introduces a small “defavorable” bias in the odds ratio estimates (and theoretically increases the p-values).

Code

tmp <- 
  col_by_week5 |> 
  bind_rows(
    col_by_week5 |> 
      dplyr::select(arm, site) |> 
      distinct() |> 
      mutate(colonized_by = TRUE),
    col_by_week5 |> 
      dplyr::select(arm, site) |> 
      distinct() |> 
      mutate(colonized_by = FALSE)
  ) |> 
  filter(! (arm %in% c("LC-106-3", "LC-106-o")))

fit <- glm(colonized_by ~ arm + site + arm:site, data = tmp, family = binomial)
fit |> summary()


Call:
glm(formula = colonized_by ~ arm + site + arm:site, family = binomial, 
    data = tmp)

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)  
(Intercept)         -2.48491    1.04083  -2.387   0.0170 *
armLC-106-7          1.87877    1.15798   1.622   0.1047  
armLC-115            2.95491    1.18673   2.490   0.0128 *
siteMGH             -0.08004    1.46978  -0.054   0.9566  
armLC-106-7:siteMGH  1.27397    1.65195   0.771   0.4406  
armLC-115:siteMGH   -0.54411    1.67176  -0.325   0.7448  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 110.619  on 83  degrees of freedom
Residual deviance:  89.848  on 78  degrees of freedom
AIC: 101.85

Number of Fisher Scoring iterations: 5

Another way is to use an actual Bayesian model that accounts for site differences, but this does not provide p-values, only posterior distributions (and confidence intervals).

Code

col_by_week5_agg <- 
  col_by_week5 |>
  group_by(site, arm) |>
  summarise(success = colonized_by |> sum(), ntrials = n())

bfit <- 
  brm(
    success | trials(ntrials) ~ arm + site + arm:site, 
    data = col_by_week5_agg, 
    family = binomial("logit")
    )


SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 3.1e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.818 seconds (Warm-up)
Chain 1:                0.901 seconds (Sampling)
Chain 1:                1.719 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 2e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.797 seconds (Warm-up)
Chain 2:                0.935 seconds (Sampling)
Chain 2:                1.732 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 3e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.728 seconds (Warm-up)
Chain 3:                1.062 seconds (Sampling)
Chain 3:                1.79 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 3e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.743 seconds (Warm-up)
Chain 4:                0.939 seconds (Sampling)
Chain 4:                1.682 seconds (Total)
Chain 4:

Code

summary(bfit)

 Family: binomial 
  Links: mu = logit 
Formula: success | trials(ntrials) ~ arm + site + arm:site 
   Data: col_by_week5_agg (Number of observations: 10) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept            -249.85    153.84  -576.84   -19.27 1.09       57      138
armLCM106M7           249.11    153.85    18.87   576.00 1.09       57      138
armLCM106M3           250.93    153.84    20.34   577.28 1.09       57      138
armLCM106Mo           252.80    153.80    22.94   579.38 1.09       57      138
armLCM115             250.46    153.85    19.76   577.88 1.09       57      138
siteMGH                10.10    273.15  -574.12   497.64 1.10       46       98
armLCM106M7:siteMGH    -8.61    273.15  -496.30   576.03 1.10       46       98
armLCM106M3:siteMGH    -9.70    273.21  -496.29   572.86 1.10       46       98
armLCM106Mo:siteMGH   471.01    331.44     3.46  1210.07 1.07       54       87
armLCM115:siteMGH     -10.92    273.17  -497.36   573.65 1.10       46       98

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Code

plot(bfit)

Code

bind_cols(col_by_week5_agg, predict(bfit))

# A tibble: 10 × 8
# Groups:   site [2]
   site  arm      success ntrials Estimate Est.Error  Q2.5 Q97.5
   <fct> <fct>      <int>   <int>    <dbl>     <dbl> <dbl> <dbl>
 1 CAP   Pl             0      11  0.00125    0.0418     0     0
 2 CAP   LC-106-7       5      15  5.00       2.48       1    10
 3 CAP   LC-106-3       8      11  7.95       2.02       4    11
 4 CAP   LC-106-o      11      12 11.0        1.30       8    12
 5 CAP   LC-115         7      11  7.00       2.15       3    11
 6 MGH   Pl             0      12  0.00225    0.0524     0     0
 7 MGH   LC-106-7       8      12  7.99       2.22       3    12
 8 MGH   LC-106-3       3       4  3.01       1.08       1     4
 9 MGH   LC-106-o       4       4  4.00       0.0354     4     4
10 MGH   LC-115         5      11  4.97       2.24       1     9

Code

plot(conditional_effects(bfit), ask = FALSE)

Code

# loo(bfit, moment_match = TRUE, reloo = TRUE)
# pp_check(bfit)

Alternatively, we can use exact fisher tests (as initially proposed by Lara Lewis, CAPRISA) for each arm x site, correcting for multiple testing, but this does not allow us to directly test for site differences.

Code

fisher_test_res <- tibble()

for (selected_site in levels(col_by_week5$site)) {
  for (selected_arm in levels(col_by_week5$arm)[-1]) {
    ft <-
      col_by_week5 |>
      filter(site == selected_site, arm %in% c("Pl", selected_arm)) |>
      mutate(arm = arm |> fct_drop()) |>
      dplyr::select(colonized_by, arm) |>
      table() |>
      fisher.test()
    fisher_test_res <-
      fisher_test_res |>
      bind_rows(
        tibble(
          site = selected_site,
          arm = selected_arm,
          odds = ft$estimate,
          pvalue = ft$p.value
        )
      )
  }
}

fisher_test_res <-
  fisher_test_res |>
  mutate(
    qval = p.adjust(pvalue, method = "BH")
  )

fisher_test_res |>
  group_by(site) |>
  gt(caption = "Fisher exact test results")

Secondary outcomes (Figures 2-4)

Kinetics of colonization

Total proportions of LBP isolates

Code

g_tot_prop_LBP_strains <- 
  colonization |> 
  as_tibble() |> 
  left_join(mae |> colData() |> as_tibble()) |>
  filter(visit_abbr != "SCR") |> 
  mutate(visit_week = visit_number - 1) |> 
  ggplot(aes(x = visit_week, y = total_prop_at, col = arm)) +
  geom_line(aes(group = pid), alpha = 0.5) +
  geom_point(size = 0.5, alpha = 0.5) +
  facet_grid(site ~ arm) + 
  scale_x_continuous("Study weeks", breaks = c(0:5,7,9,12), minor_breaks = 0:12) +
  scale_y_continuous(
    "Total proportion of LBP strains", labels = scales::percent,
    breaks = seq(0, 1, by = 0.2)
  ) +
  expand_limits(y = 0) +
  scale_color_manual(values = get_arm_colors()) +
  guides(col = "none") +
  ggtitle(make_title(mae))

Code

g_tot_prop_LBP_strains

Proportions of participants who colonized

Code

g_prop_colonized_by <- 
  colonization |> 
  as_tibble() |> 
  left_join(mae |> colData() |> as_tibble()) |>
  filter(visit_abbr != "SCR") |> 
  mutate(visit_week = visit_number - 1) |> 
  group_by(site, arm, visit_week) |>
  summarize(
    prop_who_colonized_at = mean(colonized_at), 
    # add confidence interval
    .groups = "drop"
  ) |>
  ggplot(aes(x = visit_week, y = prop_who_colonized_at, col = arm)) +
  geom_line() +
  geom_point() +
  facet_grid(site ~ .) + 
  scale_x_continuous("Study weeks", breaks = c(0:5,7,9,12), minor_breaks = 0:12) +
  scale_y_continuous(
    "Proportion of participants\nwith detected LBP strains", labels = scales::percent,
    breaks = seq(0, 1, by = 0.2)
  ) +
  expand_limits(y = 0) +
  scale_color_manual(values = get_arm_colors()) +
  ggtitle(make_title(mae)) 

g_prop_colonized_by

Code

g_prop_colonized_by

Figure 2

Code

 g_prop_colonized_by + 
  g_tot_prop_LBP_strains +
  plot_layout(guides = "collect", widths = c(1, 5)) +
  plot_annotation(tag_level = "a") &
  theme(legend.position = "bottom")

Microbiota trajectories

TODO

Code

# simplified_compositions <- 
#   mae[["ASV_all"]] |> 
#   as_tibble() |> 
#   mutate(
#     taxa = 
#       case_when(
#         str_detect(.feature, "L.crispatus[0-9]") ~ "L. crispatus (LBP strains)",
#         str_detect(.feature, "crispatus") ~ "L. crispatus (non-LBP strains)",
#         str_detect(.feature, "iners") ~ "L. iners",
#         str_detect(.feature, "Lactobacillus") ~ "other Lactobacillus spp.",
#         TRUE ~ "Non Lactobacillus spp."
#       ) |> factor()
#   ) |> 
#   group_by(taxa, .sample) |> 
#   summarize(prop = sum(ASV_prop), .groups = "drop") |> 
#   left_join(mae |> colData() |> as_tibble(), by = join_by(.sample == sample_id))

Code

# 
# g_trajectories <- 
#   simplified_compositions |> 
#   arrange_participants(by = taxa) |> 
#   filter(visit_code != 0) |> 
#   ggplot(aes(x = prop, y = pid, fill = taxa)) +
#   geom_col() +
#   facet_grid(site + arm ~ visit_code, scales = "free", space = "free") +
#   theme(panel.spacing.x = unit(0, "points")) +
#   scale_fill_manual(values = simplified_compositions$taxa |> levels() |> get_taxa_color()) +
#   scale_x_continuous(
#     "relative abundance", labels = scales::percent, 
#     breaks = seq(0, 0.5, by = 0.5), minor_breaks = seq(0, 1, by = 0.25)
#   ) +
#   scale_y_discrete("Participants", breaks = NULL) +
#   theme(
#     strip.text.y = element_text(angle = 0, hjust = 0),
#     strip.background = element_rect(color = "white")
#   )

Code

# g_trajectories

Individual LBP strains by sites

Code

LBP_strains <- 
  mae[["LBP_rel_ab"]] |> 
  as_tibble() |> 
  left_join(mae |> colData() |> as_tibble())

Code

LBP_strains_summary <- 
  LBP_strains |> 
  group_by(pid, site, arm, .feature) |> 
  summarize(max_prop = max(LBP_rel_ab), .groups = "drop")


LBP_strains_summary <- 
  LBP_strains_summary |> 
  left_join(
    mae[["mg_ksanity"]] |> 
      as_tibble() |> 
      dplyr::select(.feature, LBP, strain_origin) |> 
      distinct(),
    by = join_by(.feature)
  )

Code

# LBP_strains_summary |> 
#   filter(arm != "Pl") |>
#   mutate(strain_origin = str_c("orgin: ", strain_origin)) |>
#   ggplot(aes(x = site, y = max_prop, col = site, fill = site)) +
#   ggbeeswarm::geom_quasirandom(alpha = 0.7) +
#   geom_boxplot(alpha = 0.3, outlier.shape = NA) +
#   facet_grid(. ~ strain_origin + .feature) +
#   scale_color_manual(values = get_site_colors()) +
#   scale_fill_manual(values = get_site_colors()) +
#   scale_y_continuous(
#     str_c(
#       "Max. proportion of LBP strain\n",
#       "in active arms participants\n",
#       "at any weekly visit"),
#     labels = scales::percent, breaks = seq(0, 1, by = 0.05)
#     ) + expand_limits(y = 0) +
#   scale_x_discrete("", breaks = NULL) +
#   theme(
#     # strip.text.x = element_text(angle = 45),
#     legend.position = "bottom"
#   )

Code

LBP_strains_summary |> 
  filter(arm != "Pl") |> 
  group_by(site, .feature, strain_origin) |> 
  summarize(prop_detected = mean(max_prop >= 0.05), .groups = "drop") |> 
  mutate(origin = str_c("strain origin: ", strain_origin)) |> 
  ggplot(aes(x = .feature, y = prop_detected, fill = site)) +
  geom_col(position = "dodge", width = 0.7) +
  facet_grid(. ~ origin, scale = "free", space = "free") +
  xlab("LBP strain") +
  scale_y_continuous(
    str_c(
      "Proportion of active arms participants\n",
      "with detected LBP strain\n",
      "(rel. abundance ≥ 5% at any weekly visit)"),
    labels = scales::percent, breaks = seq(0, 1, by = 0.2)
    ) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_fill_manual(values = get_site_colors())

Introduction

Workflow

Data preparation

Loading the data and creating SummarizedExperiment objects for each assays

Clinical and survey (CRF) data

Information about the data

Loading the data

Weekly metagenomics data with LBP strains relative abundances

Checks & Quality Controls

Weekly 16S rRNA data

Checks & Quality Controls

Integrating all assays into a MultiAssayExperiment object

colData

Assays

Visit attendance

Weekly metagenomics data with LBP strains relative abundances

Weekly 16S rRNA data

Assembling the MultiAssayExperiment object

Checks, QC, and Outcomes Definitions

Checks and Quality Controls

Visualization of available data for each assay and each site x arm x participant x visit

Checking for contamination

Other checks ?

Outcomes definitions

Primary outcome definition

LBP strains relative abundances in each samples

At each visit

By each visit

Adding primary assay to MAE

Secondary outcomes definitions

16S rRNA microbiota composition, aggregated at the species level

16S rRNA microbiota composition, augmented with LBP strains relative abundances

Analyses

Demographics (Table 1)

Primary outcomes (Table 2)

Colonization by week 5

Tests / models

Secondary outcomes (Figures 2-4)

Kinetics of colonization

Total proportions of LBP isolates

Proportions of participants who colonized

Figure 2

Microbiota trajectories

Individual LBP strains by sites

Loading the data and creating `SummarizedExperiment` objects for each assays

Integrating all assays into a `MultiAssayExperiment` object

`colData`

Assembling the `MultiAssayExperiment` object

Visualization of available data for each assay and each `site x arm x participant x visit`

Adding `primary` assay to `MAE`